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To evaluate ESR spectrum at finite temperatures for specified spatial configurations of spins is 
very important issue to study quantum spin systems. Although a direct numerical estimation of the 
Kubo formula provides exact data, the application is limited to small size of the system because of 
the restriction of the computer capacity. The method of the Fourier transform of the autocorrelation 
function improved the restriction. As an extension of the method, we propose a new method for 
numerical calculation of the ESR spectrum from the time evolution of the magnetization by making 
use of the Wiener-Khinchin theorem. 

PACS numbers: 05.30.-d,75.10.jm,76.30.-v 


I. INTRODUCTION 

Quantum spin systems have attracted interests in the decades, because they exhibit various nontrivial characteris¬ 
tics. In particular, the strong quantum fluctuation and/or competition among the interaction (frustration) cause the 
singlet pair (the so-called valence bond) to behave as a unit of degree of freedom, and various novel concepts, e.g., 
the valence bond solid, resonating valence bond, and also magnon BEC, etc. have been developed. ESR is one of the 
major tools to study these properties. In particular, spatial arrangements of magnetic ions play important role for 
the property, and we need to treat the model microscopically. 

Another important topic related to the finite size effect is the effect of nonmagnetic defects in quantum spin chains. 
In spin 5 = 1/2 antiferromagnetic Heisenberg chain, the quantum fluctuation prevents the system from being ordered 
even at T = OK. A nonmagnetic defect breaks the translational symmetry and polarizes the surrounding spinsi^. 
Such a system is described by an open spin chain and some effects on susceptibility have been study on the Pd 
doped chain S^CuOg^. Due to large anisotropy no ESR signal of such defects have been report. Only magnetic 
resonance of intrinsic defects in spin-Peierls CuGeOg^ is reported since the signal from the chain drops for T < T sp . 
Recently it has been reported in organic spin chain Fabre salt, the ESR signal of an intrinsic defect in Heisenberg 
antiferromagnetic chain£. Moreover, they observed the coherence signal of the correlated defect which could have an 
important impact in the domain of quantum information2r- 10 . The knowledge of the ESR of intrinsic defect in spin 
chain is an important problem and is limited by the number of spins of the open chain one can calculate. 

To study these aspects of quantum spin systems, the direct estimation of the ESR spectrum for given Hamiltonian 
has been investigated. The most simple way is to calculate the Kubo formula directly by making use of the eigenvalues 
and eigenvectors obtained by diagonalization of the Hamiltonia n 11 1 12 . But, the application of this method is limited 
to small systems for which we can obtain all the eigenvalues and their eigenvectors. For example, for the spin system 
of N spins of S = 1/2, we need the memory of 2 2N . By making use of the symmetries of the system we may reduce 
the dimensions of the block diagonalized Hamiltonian, but still the size is limited in TV <20. 

To relieve this restriction, a method to obtain the autocorrelation function by making use of time evolution of state 
has been introduced: 13 ’ 14 The spectrum is obtained by Fourier transform of the autocorrelation function. In this 
method, The time evolution of autocorrelation function is obtained by the application of the time-evolution operator 
g-zWt ky making use of the Chebyshev iteration formula. Here we need the memory only of the order 2 N , and we 
could study double size of the case of the diagonalization method. In principle, we need to take average over the 
thermal distribution of the initial state, but thanks to the idea of the thermal stated— we do not need to take the 
complete average. By this method, the thermal property of the ESR spectrum for the single molecular magnet V 15 
which consists of 15 S = 1/2 spins has been studied^ In this method, the time evolution of the autocorrelation is 
obtained in a finite time domain 0 < t < T, and thus spectrum is suffered from the so-called Gibbs oscillation, and 
the prescription of gaussian masking is necessary to smear out the oscillation. 

In the present paper, we propose an alternate method to obtain ESR spectrum from a time evolution of a magneti¬ 
zation by making use of the Wiener-Khinchin theorem, which relates the spectrum density of magnetization dynamics 
and the Fourier transform of the autocorrelation function which gives the ESR spectrum. In quantum system, the 
definition of the dynamics of magnetization is tricky and we give a quantum version of Wiener-Khinchin relation, 
i.e, an explicit relation between Fourier transform of the autocorrelation function and the spectrum density in the 
quantum case. By making use of the relation we can obtain the ESR spectrum at finite temperatures. In this method 
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the effect of the Gibbs oscillation is significantly reduced. 

The outline of this paper is as follows. The methods previously used are explained in Sec [H] In Sec. IIII1 we 
explained the new method motivated by the Wiener-Khinchin theory. In Sec. IIVI we give the summary of the paper 
and discussion on related problems. In Appendix A, we give description of a spin system for which we demonstrate 
the methods, and in Appendix B, we explain the Gibbs oscillation. 


II. PREVIOUS METHODS 
A. Kubo formula and ESR spectrum 

The ESR spectrum is obtained by the Kubo formul a^' —. The imaginary part of the dynamical susceptibility x"( w ) 
is given by 


i r°° 

X» = 2 (! - e-^) J (M x (0)M x (t)) eq e~ iuit dt, 


and The ESR absorption spectrum is given by 


/» 




Here, we adopt usual notations: 


N 


M x (t) = e lUt M x e~ mt , M x = Sf, 

i= 1 

<-) eq = Tr[ ■ e-^]/TV[e-^]. 


(1) 

( 2 ) 

( 3 ) 

( 4 ) 


B. Numerical methods 

For numerical analysis of the ESR absorption spectrum, several methods have been developed. There are essentially 
two types of methods: (1) Exact diagonalization metho d 11 ' 12 and (2) Time-evolution of the autocorrelation function 
methodic. 


1. Exact diagonalization method 


For direct estimation of the Kubo formula, we may explicitly evaluate the formula by making use of the set of the 
eigenvalues and the eigenvectors {E n , |n )}^ =1 of the hamiltonian H where D is the dimension of Hilbert space of H: 

H\n) = E n \n) (5) 

obtained by a numerical diagonalization. The autocorrelation function (M x (0)M X (t )) eq is expressed as 


(M x (0)M x (t)) eq 


^{n\M x e lHt M x e- mt -P H \n)/Z 

n 

(n\M x \m)(m\e lEm ' t \m')(in'\M x \n)e~ lE,lt ~ l3En /Z 

n m,m' 

| (m\M x \n) feKEm-BJt-PEn /z, 

m,n 


where Z is the partition function 


( 6 ) 

( 7 ) 

( 8 ) 


z = y>-^. 


n 


(9) 
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The Fourier transform of this reads 


(M x (0)M x (i)) eq e _iwt dt 


J2\(m\M x \n)\ 2 e~^ 

m,n 


a -i (u-(E m -E r 


})t d t/Z 


J2\{m\M x \n)\ 2 e- pE ^8 (w - (. E m - E n )) /Z. 

m,n 


( 10 ) 

(ii) 


This yields the imaginary part of the dynamical susceptibility as 

1 _ p —_ 

x"M = -2- 53 lH^ x |«)| 2 e -/JBB 27 T 5 (a;-(£; r „-£; n ))/Z (12) 

m,n 

— ^ ^ (^ ^rn.n) j ( 1 * 1 ) 

m,n 

where 

£» TOin = 7r(e-^-e-^)|(m|M x |n)| 2 /^ w m , n = - E n . (14) 

We may treat only the range of u) m ,n > 0 since we are considering the absorption peak, not the emission. Note that 
X"(w) > 0 for uj > 0 . 

In this method, the spectrum is given by an ensemble of delta functions. Thus to draw the spectrum, we convert 
them in a continuous form. For example we may use bins in the to axis, or we replace the delta function by a gauss 
distribution with some variance which represents a finite resolution. 

This method is exact, but we need to obtain all the eigenvalues and their eigenstates. Therefore, we need to store 
the matrix of the size D, which is 2 N for systems of N spins with S = 1/2. This requires D 2 in the memory. If 
the system has symmetry, we may reduce the size. For example if the system conserves M z , then D is reduced to 
nCn/ 2 +m z ■ Moreover for ESR only the uniform mode is relevant and thus only the fully symmetrized states are 
necessary, which also reduces D. However, the limitation of the memory prevents us from calculation more than 
N = 20 for systems of S = 1/2. In Fig. |T] an example of spectrum obtained by exact diagonalization method for an 
antiferromagnetic Heisenberg chain with N = 14 with the static field H = 5K at the temperature T = 100K. The 
notation of the model is explained in Appendix A. The data is given as ensemble of the delta peaks (Eq. (fT3l) l. we 
made a histogram with small bins to show the spectrum with a mesh of the frequency Aw = 0.00005. 


£ 



5.001 5.002 5.003 5.004 5.005 


FIG. 1: The spectrum of an antiferromagnetic Heisenberg chain obtained by exact diagonalization method: N = 14, T = 100K, 
H = 5K, and Au = 0.00005. 
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2. Time-evolution of the autocorrelation function method 


As we find in the Kubo formula(Eq. ([]}), the spectrum is given by the autocorrelation function. Thus, if we have 
the autocorrelation function, the spectrum is given by its Fourier transform. 

First, we derive the explicit expression to calculate in the time-domain method. The integral is divided to two 
parts. 



(M x M x (t)) eq e~ iwt dt 


(M x M x (t)) eq e~ lut dt - 


(M x M x (t)) eq e~ iut dt. 


The second integration term can be reduced to the following form: 



{M x M x (t)) eq e~ iut dt 


pOO 

/ (. M x M x (-t)) eq e +iut dt 

Jo 



((M x M x (t)) eq e~ iujt )* dt, 


(15) 


(16) 


where the first equal sign follows from the change of variables t —> — t and the second one the relation 
(( M x M x (t )) eq )* = (M x M x (—t)) eq (* denote complex conjugate). Therefore we have 


pOO 

(M x M x (t)) eq e~ iut dt = / {M x M x (t)) eq e~ iuit + ({M x M x (t)) eq e~ iut Y 

3 Jo L 


dt 


= 2 Re 


(M x M x (t)) eq e~ iut dt 


(17) 

(18) 


where Re[-] denotes the real part. 

The autocorrelation function ( M x M x (t)) eq is given by 


(M x M x (t)) eq = 


TrM x e mt M x e- im e-P' H 

Tre"^ 


(19) 


For a small system whose eigenstates can be obtained, we can calculate the trace numerically exactly. However, for 
larger systems we can obtain it by making use of what we call random vectors^* —. In this method we prepare the 
so-called Boltzmann-weighted random vectors or typical state (see below)^r— 1$^), and perform time evolution of 
the state by applying which can be done by the Chebyshev method. In this method, we need memory only of 

the order D not D 2 . Thus we can calculate up to larger number of spins. 

Let us explain the method of the typical state briefly. Let {|n)}^L 1 be an arbitrary set of complete orthonormal 
states of Hilbert space. Using the complex-valued random variables {Cnj^Li, we define the random vector as 


D 

b) = J2tn\n). 

n —1 


The coefficients {£n}n satisfy the following conditions: 


E[&J = 0, 

E[£m£n] = 0, 

E [£Un] = Smn- 


( 20 ) 


( 21 ) 

( 22 ) 

(23) 


In this paper, we draw them from the 2D dimensional spherical surface: X)«=i l£™| 2 = -D- This construction is 
independent of the choice of the set of bases {|n)} n . The thermal state (or typical state) in the system Hamiltonian 
TL at an inverse temperature /3 is defined by 

|^>=e-3/”». (24) 


Using this state, we can calculate the expectation value of an observable X in equilibrium by 

Tr[Xe~^] EK^IXI^)] 

Weq Tr[e-W] E[(^|^)] ' 


(25) 


In the program, we take the average E[-] with respect to the same samples {| < 1>^)} in the dominator and in the 
numerator. In fact, it is known that the convergence of E[-] becomes faster with the increase of the dimension D of 
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the Hilbert space. If we evaluate this quantitatively by using finite samples (j = 1, • • ■ S) of the states the deviation 
is given 


( |Ti'[A'] - § £?=i(^|A1^>| 2 ^ \ ^ 1 DTr[X*X] - \Tr[X]\ 2 

{ \Tr[X}\ 2 - a )-aS(D+ 1) \Tr[X)\* 


(26) 


So it would be enough to prepare only a few typical states as samples for a sufficiently large system. 
The autocorrelation function ( M x M x (t)) eq is given by 


(M x M x (t)) e q = 


E[(® fl \M x e mt M x e- i ' Ht \‘f>f,)] 


The essence of this method is how to deal with e and e Mt . 

Here, we review this problem introducing the expansion with the Chebyshev polynomial as follows: 


(27) 


|$/3) 


e /3A/2 e rHsc \v) 


-0X/2 


Io(-t)T 0 (H sc ) + 2 2 I k (-T)T k (H sc ) 


k= 1 


l$/3>, 


where 


U = A\n sc + A, AA 


A = 


f'max + -Err 


(28) 

(29) 


(30) 


and E? max and E m - ln are the largest and smallest eigenvalues of TL. The infinite sum should be truncated after & max 
terms, which is chosen such that Efc max (—r) is sufficiently small. //.(•) denotes the modified Bessel function of order k 
and T/ c (■) the Chebyshev polynomial of order k. The Chebyshev polynomials satisfy the following recurrence relation: 


T k+1 (H sc )=2H sc T k (H sc )-T k _ 1 (H sc ), T 0 {H SC ) = 1, T 1 (H BO )=H BC . (31) 

Using this expansion and recurrence relation, we can obtain |<l> j a) only multiplying a vector by "H sc repeatedly, and 
summing them up, without storing any large matrices such as the Hamiltonian "H. We can also treat the time evolution 
operator e~ lHt similarly. 

In this way, we obtain the time-series data of the vectors: 

M x e - iw \<f>p), and e~ mt M x \$p) (32) 

and the autocorrelation function f(t) = (<h p\M x M x (i)|$ p). 

Next, we apply discrete Fourier transform(DFT) on them. Specifically, using the set of discrete values {f k } k Zo = 
{/(^■)}fc=o f° r gi ven T and n, we obtain the DFT result as a set of discrete values 



where only the first half of { ui k } k are significant since the latter half correspond to the negative frequencies, and the 
maximum frequency is less than . It should be noted that here the range of the time integral is finite but not oo as 
in the definition of original Fourier transform. This leads to what is called the Gibbs oscillation problemi^ (Appendix 
B). The Gibbs oscillation with the negative values of the spectrum occurs due to the discontinuity at the ends of the 
interval of integration. To avoid this phenomenon, we may apply a suitable window function to DFT—. Here, we use 
a Gaussian window: 


[ /(t)e- iwt e-K a T) 2 dt = 2Re [ dt 

J-T JO 


(34) 


The Gaussian window function smears the Gibbs oscillation and gives a Gaussian-like profile with correct amplitude. 
As depicted in Fig. [2] the Gibbs oscillation disappears and the result is consistent with the one derived with the exact 
diagonalization method with the resolution of 0(a). 













FIG. 2: The spectrum obtained by exact diagonalization method and by the time-evolution method : N = 6, temperature 
T = 500K, H = 5 K, sampling time dt = 0.5, the number of data n = 100000, and the mesh of the frequency A ui = 0.00005 (red 
line) and 27r/50000 (blue line), (a) the spectrum obtained by exact diagonalization method (red line) and by time-evolution 
method without window functions (blue line). You can see the negativeness of the spectrum caused by the Gibbs oscillation, 
(b) The smeared Gibbs oscillation owing to the Gaussian window function with a = 6 which means that the variance in the 
frequency domain is ^ = 0.00012. Compared with the spectrum obtained by exact diagonalization with the delta peaks replaced 
by the Gaussians with variance = (0.00015) 2 (red line), the time-evolution spectrum (blue line) is consistent as expected. 


III. WIENER-KHINCHIN METHOD 

Now, we propose a new method making use of the Wiener-Khinchin theorem. In this method, we use the time- 
evolution of the state as in the previous subsection, and thus by this method we can study large systems as well. But 
here, instead of the autocorrelation function, we study the dynamics of magnetization itself, i.e., ( M x (t )) from an 
initial state 1$^), and use the Wiener-Khinchin relation, that is, the relation between the Fourier transform of the 
fluctuation in time of a quantity A(t), and the spectrum of the autocorrelation function (X(0)X(t)). 


A. Wiener-Khinchin theorem 


First, we briefly review classical Wiener-Khinchin theorem. For a time-evolution of a quantity X(t), we define the 
autocorrelation function R(t) and the spectral density S(uj) as 

1 f T / 2 

R(t) = (X(0)X(t)) = lim - / X(r)X(t + r)dr, (35) 

T~T J _ T/2 


and 


l* T M | 2 


! r T / 2 

S(“) = Jim r' /l , X t (uj)= X{t)i 

T-*oo 1 J — T/2 


—i cot 


dt, 


(36) 


respectively. The Wiener-Khinchin theorem tells us that the Fourier transform of the autocorrelation function equals 
to the spectral density: 


G{u)= / R{t)e~ iuit dt = S{u). 


(37) 


Note that here we assumed that process X(t) is stationary. 

The autocorrelation function R(t) is defined as the autocorrelation function along the time-evolution. We also 
assume that the average is the same as that in the ensemble average in the equilibrium state if the process X(t) is 
ergodic. 
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B. Dynamics of the magnetization 

Now we go back to the original model in question. We need the autocorrelation function R(t) = (M x (0)M x (t)) eq to 
obtain the spectral density S(u>) which is directly related to the ESR spectrum. Here it should be noted that M x (t) 
is an operator and the definition of ( M x (t )} is tricky. Definitely, ( M x (t)) eq is time-independent. In the following, we 
propose a numerical method motivated by the Wiener-Khinchin theorem, but it should be understood as a numerical 
algorithm to obtain S(uj)M 

First we prepare a typical initial state for the canonical ensemble, 1$^} = e~^ H ^ 2 \v). The expectation value with 
the state gives the equilibrium state at the given temperature. But for an sample of the typical state, ( M x (t )} is not 
necessarily zero. Then we calculate the time evolution of the selected state and obtain the value of ( M x (t )} at each 
time. From this time evolution of ( M x {t )}, we obtain the spectral density Eo. flMl) . 

At last, by averaging the spectral densities with respect to several initial states, we obtain the spectral density 
for this process. Here we again note that the expectation value E [(M x (t))] in the equilibrium is zero. Thus, the 
expectation value of the spectral density of its Fourier transform 

fT /2 

Mj{uj)= / {M x (t))e~ iut dt (38) 

JT/2 


is also zero. 

However, the expectation value of the spectral density E[|Mj(w)| 2 ] is not zero, and this provides information for 
the Fourier transform of the autocorrelation function G(u>) for the ESR spectrum of the system. We will show below 
the explicit relation of this spectral density and G(uj). 


C. Relation between the spectrum density and G(u) 

The time evolution of the system from a given initial state |$^) provided with the random coefficients {£„} is 
assumed to be given by the quantum mechanical evolution e~ lW . 

Thus, we introduce a quantity M|(t) by 

Mf(i) = (^|AP(t)|^) = Y,^ne-^ Em+E ^/ 2 e i ^- E ^ t (m\M x \n), (39) 

m,n 

where | n) and E n are the eigenvector and its eigenenergy of the system Hamiltonian R. 

The Fourier transform Mj (w) is expressed as 

Mj( w ) = Y,Cmtne- f){Em+En)/2 2TrS T (to - (E m - E n )) (m\M x \n), (40) 

m,n 

where 

-i pT /2 oin 

S T (to) = e~ iut dt = -3- -1=^ 6(w). (41) 

27r J-T/2 


Therefore 

|MjM| 2 =]T ^ Crn^rn'Cn'e- 0 iEm+En)/ 2 e -^' +E ^ /2 

m,n m' ,n' 

x dir 2 S T (w - (E m - E n )) 5 t (w - (E m , - E n ,)) (m\M x \n)(n , \M x \m'). (42) 

Here we take the average over the random variables {£„} n . Using the following formula^: 


E[C£n£m'&] = — ^ ,n' (1 ) H - (1 &m,n )}+ dTT' 5 m,n' 3n,m' ^m,n ? 


and (n\M x \n)=Q, we have 


E[|MjM| 2 ] = 5^E[|^| 2 |^| 2 ]e- /3(£7 " 1+i57 - ) 47r 2 (<5 T (w - (E m - E n ))f \(m\M x \n)\ 2 

m,n 


( 43 ) 






(44) 


= p^jTe-P“J2 e ~ 2flEn2n6T ( w - - E ^))\(m\M x \n)\ 2 , 

m,n 


where we used the relation 


(S t (uj - (E m - E n ))) 2 k, ^-5 t (uj - (E m - E n )), 


and the hnal result for the spectral density reads 

e[|mJmi 2 ]/z| 


T,g(co) = lim 
r T-> oo 

D 

D + l 


T 


e -0" J2 e~ 2 P En 2nS (w - (E m - E n )) \(m\M*\n)\ 2 /Z% 


(45) 


(46) 

(47) 


where 


Zp = E[(^|^)]. (48) 

Finally by comparing this result with the Fourier transform of the autocorrelation function 

Gp{u>) = 53e-^»27r5(w - (E m - E n ))\{m\M x \n)\ 2 /Z p , (49) 

m,n 

we obtain a Wiener-Khinchin-like relation 

Gp(u) = ^i^e^S^ /2 ( w ). (50) 

L) Zip 

Thus, the imaginary part of the dynamical susceptibility of our interest is given by 

= ~^~~z^ sinh (Jy) ^p/ 2 ^' 

Here it should be noted that we need to calculate the quantities of /3/2 (not 0) to obtain the ESR spectrum of 0 
In the present method, we obtained the approximate spectrum from a finite time domain. The finiteness of the 
domain causes the Gibbs oscillation in this case as well. However, because we take the square of Mj(w) fEn. (|42]l ). 
and thus the negative value would not appear. 

Here we note that the total amplitude of the spectrum, i.e., the integration over the spectrum does not change 
because of the relation: 


1 

T 



^sinwT/2^ 



(52) 


and we can obtain the correct spectrum in the limit T —> oo in the present method. 

As to the effect of the Gibbs oscillation, in the present method it gives the width of the spectrum due to the 
finiteness of the observation. In this sense we regard that the width comes from the way of observation and is natural, 
in contrast to the case of the method presented in the section III B 21 where we introduced a gaussian window with an 
artificial width to smear the Gibbs oscillation which is the same order of the natural width. The width of the gaussian 
window is taken to be in the same order of the width of the Gibbs oscillation. 

We show a comparison of the spectrum obtained by the exact diagonalization method and by the Wiener-Khinchin 
method. 
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FIG. 3: Absorption spectrum obtained by the exact diagonalization method and by the Wiener-Khinchin method: N = 6 (100 
samples), T = 200K, H = 5K, and the mesh of the frequency = 0.000063 (red line) and 27r/100000 (blue line). 

We also show the intensity 


f‘00 

r= r{u)duj (53) 

Jo 

which is the integral over the spectrum I x (uj ) obtained by the exact diagonalization method and by the Wiener- 
Khinchin method as a function of temperature. The data of WK method are obtained by averaging 100 samples. 
Here we find almost perfect agreement. 



T( K) 


FIG. 4: Comparison of intensities of the absorption as a function of the temperature: N = 6 ,H = 5K, the number of samples 
S = 100. 
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D. Typicality of the fluctuation 


If Tr[X] ^ 0, the estimation Ea. (l26l) works. But, in WK method Tr[Xe _,3W ] = 0, and thus we cannot use this bound 
for the estimation of the error for the quantity Mf(f) (Eq. (l39l) l or Mf(w) (Eq. (1401) 1. Thus to obtain E[|Mj(w)| 2 ] 
(Eq. (14411 1 from numerical data, we need to estimate the variance of the quantity: 


E 


(|Mj( W )| 2 -E[|MjM| 2 ]) 2 


= E[|MjH| 4 ]-E[|MjM| 2 ] 2 


(54) 


If we assume the non-degeneracy of energy gaps that if E m — E n = E m / — E n / then (m = m! and n = n') or (m = n 
and m' = n’), and approximate 5 T (x ) = S(x ), then we have 

E ->• E [(C( m ) 2 K(n) 2 ] - 4, (55) 

and 

E[|Mj(o;)| 4 ]-E[|Mj(a;)| 2 ] 2 ^3E[|Mj(a;)| 2 ] 2 . (56) 

Thus, the distribution of the quantity | Mj (w)| 2 has a standard deviation 

v^E[|MjM| 2 ], (57) 

which does not depend on D and of the order of the value of average. This means that the distribution function 
converges to a fixed form. Thus, we can estimate E[|Mj(w) | 2 ] by some ensemble average over the sample independently 
of D. 

In the case of the autocorrelation function Tr[A'e _,3W ] ^ 0, and we expect that the distribution of the obtained values 
converges to the expectation value with the variance of the order 1/D. In this sense, the method of the autocorrelation 
has an advantage. But the WK method can also give a value with finite sampling because the variance of which does 
not depend on D, and the distribution shows a kind of typicality 2 ^ We may call it the typicality of the distribution 
for zero-mean quantity. 


IV. SUMMARY AND DISCUSSION 

We have proposed a time-domain method to obtain ESR spectrum by making use of Wiener-Khinchin relation, 
in which effects of the Gibbs oscillation due to finite observation period is practically suppressed. This method is 
used the thermal typical state for the variable whose expectation value is zero in the thermal average, and thus the 
temperature of the thermal state and that of the spectrum differ by the factor 2. 

Here we presented one of the way for time-domain method motivated by the Wiener-Khinchin theorem, but there 
are many other ways to obtain the spectrum density by making use of relation related to Wiener-Khinchin theorem. 
For example, if we calculate E[(<b j g|M a '| < h / 3 )(<b / 3 |M :r (f)| , h ( 3 )], we can obtain a similar expression. The relation among 
them will be interesting problem. 

It should be noted that the present method does not directly relate to the experimental situation. In this paper, we 
proposed the present method by making use of mathematical relations discussed in Sec 1III1 to obtain the spectrum. 
But, dynamics from the thermal typical state is an interesting problem to study, which will be studied elsewhere. 
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Appendix A: Model for Application 


We demonstrate the methods applying to the one-dimensional spin-1/2 XXZ model with the dipole-dipole interac¬ 
tion under a static magnetic field H and an oscillating field A(f). We obtain the response to the AC field A(t) along 
the x axis. The hamiltonian of the system is given by 


where 


Htot = n + \(t) 

= Ho + n' + Hn + Hz + \(t), 


N -1 


n 0 = J Y, S > ■ s i+ 1: 


i=1 
N—l 


w = 

2—1 

Hb = Do 


Si-S : 3 (Si ■ r. ij )(S j ■ r^) 


r S 

(i,j) V ij 

N 

Hz = -9HbH Sf, 
2=1 
N 

A (t) = ApCOS 


r- 

ij 


(Al) 

(A2) 

(A3) 

(A4) 

(A5) 

(A6) 

(A7) 


Hereafter we put g^B = 1- We adopt Kelvin as the unit of enegy, and A = 0.00001. For the demonstration we put 

Do = 0 . 


Appendix B: Gibbs oscillation and window function 


Here we discuss the properties of Gibbs oscillation which inevitably comes from the finite observation time. This 
is a general problem encountered in performing the Fourier transform of time-series data in finite time domain^. In 
the method of autocorrelation function this effect causes undesirable negative values of the spectrum, which has been 
managed to be smeared by the use of a window function. 

First let us have a brief review on the Gibbs oscillation. Let {fk}kL -oo be the time sequence of data such as the 
autocorrelation function. The spectrum of process is given by 


F{v) 


fk 


00 / f‘00 \ 

A f* e ~ iUkA (- f(t)e- iut dt), 

k =-oo ^ J 


/*7r/A 

- / F(u;)e 


iujkA 


2n 


doj, 


’ — 7r/A 


(Bl) 

(B2) 


where A is the interval of sampling. This transformation is often called the discrete time Fourier transform (DTFT). 
In reality, however, since we cannot prepare an infinite number of data, we have to terminate the sum at a finite 
number. Thus the transformation is modified as 


n—l 

F(w) = A ^2 f k e 

k=—n 


—iuik.A 


oo 


A ^ fkh k e 11 


(B3) 


h k — 


(B4) 


where we used a symmetry of f k 


1 

0 


k = 0,±l,±2...,±(n- 1), —n 
otherwise. 
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This formula is expressed in the form of convolution: 


F M 


1 

27r 



F(u')H(u 




1 

27T 





(B5) 


where -ff(w) is the DFTT of {/ifc}fc 


n—1 

H(u) = A ^ e ~iukA 

k=—n 


= Ae 1 ' 


: sin(omA) 


sin(^) 


(B6) 


Above formula means that the spectrum we need is deformed by the DTFT of the rectangular window H(uj). In 
Fig. 0a), H(uj) = A sin(uA/ 2 ) depicted. This oscillation gives the negative peak of the spectrum F(lo). (Note that 
the frequency ui is also discretized in the finite domain.) 

In order to avoid this apparent negative peaks, a window function method unusually with a Gaussian window has 
been introduced. 


9k = h k e~^ a ™ )2 , k = 0, ±1, ±2,... (B7) 

In Fig. [5](b), the spectrum of |G(u;)| is depicted. 

The parameter a determines the artificial resolution of the spectrum. Within this resolution, the Gibbs oscillation 
is smeared out and, so we can reproduce the spectrum as discussed in section III B 21 



- (a) 

- (b) 


FIG. 5: Square root of power spectrum of window function : n = 10000, A = 0.5, and a = 6. (a) rectangle window, (b) 
gaussain window. 

Here, we take the value of A which is the interval of observation time arbitrary. But, in the case the eigenenergy 
is confined in a finite range E m [ n < Ei < E max , that is the situation so-called Bond-limited function, the Nyquist- 
Shannon theorem tells us that A should be smaller than n/E max . Thus, it is most efficient to take A in the Chebyshev 
procedure to be this value. 

The width of the Gauss window a is known to be taken as 

e -K“S) 2 = £ ^a 2 = -21n£, (B8) 

where e is a number of the order of the smallest number of the computer resolution, say e = 10 -12 . This is consistent 
with the above choice of a = 6. The period of the Gibbs oscillation is of the same order of the width of the gaussian 
window. 

In contrast to the method using the autocorrelation function, in the WK method the spectrum is squared, and 
the spectrum is given by |T(w)| 2 and thus the spectrum is positive. Although the effect of the Gibbs oscillation 
exists, it gives natural width due to the finite time window and do not harms the spectrum much. In Fig. [6l we show 
comparison of H(ui) and \H{lo)\ 2 
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- (a) 

- (b) 


FIG. 6: Comparison of windows : n = 10000, A = 0.5 : (a) (b . 
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